Labwork 3 : Estimation for PALM microscopy

Gabriel Gostiaux

Contents

1. Location estimation
2. Nuisance parameters
3. Influence of the imaging system characteristics

---

Due to light diffraction, the resolution of conventional microscopic images is limited by the wavelength. Nevertheless various methods have been developed in order to "beat" the diffraction limit, especially in fluorescence microscopy. One of these methods is the photoactivated localization microscopy (PALM). This technique consist in isolating the observed fluorescent emitters and fitting their PSF. This way one can determine the location of the isolated emitter with a precision better than the Rayleigh criterion, achieving "super resolution". In this labwork, it is proposed to estimate the location of a given emitter.

For simplicity reasons and without any loss of generality, we consider the location from a 1D point of view. Let $x$ be the spatial coordinate, $\theta$ the position of an emitter. It is then assumed that the signal given by this single emitter on the image plan can be written as $s(x, \theta)=a r(x, \theta)+b$ where $r$ is a Gaussian function centered on $\theta$ so that

$$ r(x, \theta)=\frac{1}{\sqrt{2 \pi} \sigma_{r}} \exp \left[-\frac{(x-\theta)^{2}}{2 \sigma_{r}^{2}}\right] $$

where $\sigma_{r}$ is expressed as a function of the full width at half maximum (FWHM) $w$ of the signal on the sensor. In this labwork, $b$ is assumed to be an additive normal white noise of standard deviation $\sigma_{b}$.

P1 We define $r_{i}$ the integration of $r(x, \theta)$ over the pixel $i$ with $i \in[0, N-1]$

$$ r_{i}=\int_{i \Delta x}^{(i+1) \Delta x} r(x, \theta) \mathrm{d} x $$

where $\Delta x$ is the discretization step. Give the expression of $r_{i}$ as a function of the parameters $\theta$ and $w$. We remind that

$$ \begin{equation*} \frac{2}{\sqrt{\pi}} \int_{0}^{z} e^{-u^{2}} \mathrm{~d} u=\operatorname{erf}(z) \tag{3.1} \end{equation*} $$

Let's compute $r_i$ with a integration by substitution :

...

---

1 Location estimation

We assume in this section that $a$ is known and $a=1$.

P2 Give the expressions of the CRLB and the Maximum Likelihood (ML) estimator for the estimation of $\theta$.

...

$\rightsquigarrow \quad$ Assuming $N=100, w=2 \mu m, \Delta x=2 \mu m$, plot the the signal $r\left(x, \theta_{0}\right)$ versus $x$ for a few values of $\theta$ (for example $\theta=\theta_{0} \in\{46.0, 47.2, 48.4, 49.6, 50.8, 52.0\}$ (in $\mu m$ ).

$\rightsquigarrow \quad$ On the corresponding figures, plot simultaneously the signal $r\left(x, \theta_{0}\right)$ and $r_{i}$ for a few values of $\theta$. You will plot each curve on a different figure using for example subplot and bar Matlab functions. For these figures, you can reduce the $x$-interval to $[40 \mu \mathrm{~m}, 60 \mu \mathrm{~m}]$.

$ i \in [0,99] \quad x \in [0,198] $

import numpy as np
from math import pi, sqrt
import plotly.graph_objects as go
import scipy as sp

import plotly.io as pio
pio.renderers.default = 'notebook'  # For classic Jupyter Notebook

N = 100
w = 2
dx = 2
sigma_r = w /2 /sqrt(2*np.log(2))

L_theta_0 = [46.0, 47.2, 48.4, 49.6, 50.8, 52.0]

def r(x, theta_0):
    """Compute the localization of a point according to a gaussian function of mean the true position.

    Args:
        x (list): position
        theta_0 (scalar): true position

    Returns:
        list: localization
    """
    x = np.asarray(x)

    res = 1 /sqrt(2*pi) /sigma_r * np.exp( - np.power(x - theta_0, 2) /2 /sigma_r**2 )

    return res


def r_i(x, theta):
    
    x = np.asarray(x)
    res = 1 /2 * ( sp.special.erf((np.add(x, np.ones(x.shape)) * dx - theta) /(sqrt(2) * sigma_r)) - sp.special.erf(((x) * dx - theta) /(sqrt(2) * sigma_r)) )

    return res


def run():

    L_x = [i for i in range(N)]

    # print(f"Start : {L_x[0]} - Stop {L_x[-1]} - len : {len(L_x)}")

    L_r = []
    L_R_i = []

    # Create the plot
    fig = go.Figure()

    # Add error bars
    for t in L_theta_0:
        R = r(L_x, t)
        R_i = r_i(L_x, t)
        L_r.append(R)
        L_R_i.append(R_i)

        fig.add_trace(go.Scatter(
            x=L_x,
            y=R,
            mode='markers+lines',
            name=f'Localization for theta_0 = {t}',
            #error_y=dict(type='data', array=std_devs, visible=True),
            marker=dict(size=4)
        ))

        fig.add_trace(go.Scatter(
            x=L_x,
            y=R_i,
            mode='markers+lines',
            name=f'Discret localization for theta_0 = {t}',
            #error_y=dict(type='data', array=std_devs, visible=True),
            marker=dict(size=4)
        ))

    # Customize layout
    fig.update_layout(
        title='Localization for different true positions theta_0 ',
        xaxis_title='position',
        yaxis_title='localization',
        #yaxis=dict(title='LLHs', range=[0, max(LLHs) + 1]),
        #xaxis=dict(title='Decay Rate (alpha)', tickvals=decay_rate),
        showlegend=True
    )

    # Show the plot
    fig.show()


if __name__ == "__main__":
    run()

$\rightsquigarrow \quad$ On another figure, plot the CRLB for the estimation of $\theta$ as a function of $\theta$. You can choose $\sigma_{b}=0.01$.

sigma_b = 0.01

def noisy_R(x, theta):

    n_R = r(x, theta) + np.random.gaussian(0, sigma_b, shape=x.shape)

    return n_R

def CRLB(x, theta):
    
    a = 1

    res = sigma_b**2 /sigma_r**2 /a /np.sum(np.exp( - ((np.add(x, np.ones(x.shape)) * dx - theta) /(sqrt(2) * sigma_r))**2 ) - np.exp( - ((np.add(x, np.ones(x.shape)) * dx - theta) /(sqrt(2) * sigma_r))**2 ))

    return res

Q1 Comment the shape of the CRLB curve, its dependence to $\theta$. Considering the behavior of $r_{i}$ with $\theta$, propose an interpretation about the shape of the CRLB.

$\rightsquigarrow \quad$ Design a Matlab script or function to generate a realization of the noisy signal.

$\rightsquigarrow \quad$ Design a Matlab script or function which finds the maximum likelihood estimator from a realization. You can use the Matlab functions fzero or fminsearch.

$\rightsquigarrow \quad$ Estimate the statistical mean and variance of the ML estimator with a Monte-Carlo simulation.

Q2 Can we say that the estimator is efficient? Justify your answer.

Useful commands/tips:

• $\operatorname{bar}(\mathrm{X}, \mathrm{Y})$ : draws the columns of the M-by-N matrix Y as M groups of N vertical bars. The vector X must not have duplicate values.
• $\operatorname{subplot}(\mathrm{n}, \mathrm{m}, \mathrm{k})$ : divides the current figure window in $n \times m$ subfigure and considers the $k^{\text {th }}$ subfigure. The subfigures are numbered from left to right and from the top to the bottom.
• erf: is the error function defined by equation (3.1)
• hold on: holds the current plot and all axis properties so that subsequent graphing commands add to the existing graph.
• fzero: tries to find a zero of a function. The syntax to use is : fzero(@(x) fonc(x,a_1, a_2), X0) to find the zero of the function fonc( $\alpha$, a_1 , a_2) according to $\alpha$. The a_1, a_2 parameter stay constant during optimization. XO is the initial value.
• fminsearch: attempts to find a local minimizer of a function. The syntax to use is : fminsearch(@(x) fonc(x, a_1, a_2), X0);

2 Nuisance parameters

We now assume that $a$ is unknown, since it is of no interest for the location of the emitter it is considered as a nuisance parameter.

P3 Propose a method to deal with this kind of parameters.

We can select it as a parameter and first etimate it to lower the CRLB of $\theta$.

P4 Give the expressions of the CRLB and the ML estimator for the estimation of $\theta$ in this case.

$\rightsquigarrow \quad$ Plot the CRLB of the estimation of $\theta$ as a function of $\theta$.

Q3 Compare it with the case where $a$ was known: what can you tell? Was this result expected and why?

$\rightsquigarrow \quad$ Design a Matlab script or function which finds the maximum likelihood estimator from a realization in this case.

$\rightsquigarrow \quad$ Estimate the statistical mean and variance of the ML estimator with a Monte-Carlo simulation.

Q4 What are your comments?

3 Influence of the imaging system characteristics

We still assume that $a$ is unknown. The width $w$ can be seen as a parameter that depends on the imaging system characteristics, as it can be tuned by changing the magnification of the system. We propose here to find the value of $w$ that will lead to the best location estimation.

$\rightsquigarrow \quad$ Plot the CRLB calculated in section 2 for various values of $w$ (for example $w \in\left\{\begin{array}{llll}1 & 2 & 3 & 4\end{array}\right\}$ ).

Q5 Comment the shape of the curves. How does the CRLB evolve according to $w$ ?

$\rightsquigarrow \quad$ Plot the evolution of maximum of the CRLB as a function of $w$.

Q6 Comment the curve. How would you choose the magnification of your imaging system to have the best location precision? What is the expected resolution in this case?